void ASTroCAT_chemostat_expts() {
  Environment env;
/*  
  env = new ASTroCAT_Chemostat(2);
  env.resetParams();
  env.getParam("do_PARcycle").current = 1;
  ac_chemostat_interannual_run(env, "interann_np2");

  env = new ASTroCAT_Chemostat(4);
  env.resetParams();
  env.getParam("do_PARcycle").current = 1;
  ac_chemostat_interannual_run(env, "interann_np4");

  env = new ASTroCAT_Chemostat(8);
  env.resetParams();
  env.getParam("do_PARcycle").current = 1;
  ac_chemostat_interannual_run(env, "interann_np8");

  env = new ASTroCAT_Chemostat(16);
  env.resetParams();
  env.getParam("do_PARcycle").current = 1;
  ac_chemostat_interannual_run(env, "interann_np16");

  env = new ASTroCAT_Chemostat(32);
  env.resetParams();
  env.getParam("do_PARcycle").current = 1;
  ac_chemostat_interannual_run(env, "interann_np32");
*/

  env = new ASTroCAT_Chemostat(64);
  env.resetParams();
  env.getParam("do_PARcycle").current = 1;
  ac_chemostat_interannual_run(env, "interann_np64");

  println("done");
}




void ac_chemostat_run(Environment env, float intLength, String name) {
  float saveTimestep = 1;
  env.updatePrivateParams();
  env.rewind();
  env.startArchive(saveTimestep);
  env.archive.saveFluxes = true;
  env.archive.exclude("grazP*");
  env.archive.basename = name;
  println("running " + name);
  float lastReport = 0;
  while (env.time.current[0] < intLength) {
    env.integrate(env.archive.dt);
    if ((env.time.current[0] - lastReport) > intLength / 20.) {
      lastReport = env.time.current[0];
      println("    t = " + env.time.current[0]);
    }
  }
  env.saveArchive();
}


void ac_chemostat_interannual_run(Environment env, String name) {
  float saveTimestep = 1;
  int nYears = 50;
  float q_mean = 0.03;
  float q_std = 0.005;
  randomSeed(1);
  env.updatePrivateParams();
  env.rewind();
  env.startArchive(saveTimestep);
  env.archive.saveFluxes = true;
  env.archive.exclude("grazP*");
  env.archive.basename = name;
  println("running " + name);
  for (int i=0; i<nYears; i++) {
    println("    year " + i);    
    env.getParam("q").current = q_mean + q_std * randn();
    float intLength = 365 * (i+1);
    while (env.time.current[0] < intLength) {
      env.integrate(env.archive.dt);
    }
  }
  env.saveArchive();
}







void ASTroCAT_upwelling_expts() {
  int NP = 2;
  float Ntot = 40;
  Environment env = new ASTroCAT_Upwelling(4, 8, 4e4, 8e4, NP, Ntot);
  env.resetParams();
  String basename = "ac_upw_np" + NP;
  
  ac_upw_run(env, basename);
  env.getParam("zeta").current = 4;
  ac_upw_run(env, basename + "_zeta4");
}



void ac_upw_run(Environment env, String name) {
  float intLength = 180;
  float saveTimestep = 1;
  env.updatePrivateParams();
  env.rewind();
  env.startArchive(saveTimestep);
  env.archive.basename = name;
  println("running " + name);
  float lastReport = 0;
  while (env.time.current[0] < intLength) {
    env.integrate(env.archive.dt);
    if ((env.time.current[0] - lastReport) > intLength / 20.) {
      lastReport = env.time.current[0];
      println("    t = " + env.time.current[0]);
    }
  }
  println("    saving");
  env.saveArchive();
  println("    saved");
}
